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Abstract 



For a large class of processes with an absorbing state, statistical 
properties of the surviving sample attain time-independent values in 
the quasi- stationary (QS) regime. We propose a practical simulation 
method for studying quasi-stationary properties, based on the equa- 
tion of motion governing the QS distribution. The method is tested in 
applications to the contact process. At the critical point, our method 
is about an order of magnitude more efficient than conventional simu- 
lation. 
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Stochastic processes with an absorbing state arise frequently in statistical physics 
[1,2], epidemiology [3], and related fields. In autocatalytic processes (in lasers, chem- 
ical reactions, or surface catalysis, for example), and population models, a popula- 
tion of n > individuals goes extinct when the absorbing state, n = 0, is reached. 
Phase transitions to an absorbing state in spatially extended systems, exemplified 
by the contact process [4,5], are currently of great interest in connection with self- 
organized criticality [6], the transition to turbulence [7], and issues of universality in 
nonequilibrium critical phenomena [8-10]. Such models are frequently studied using 
deterministic mean-field equations, Monte Carlo simulation, and renormalization 
group analyses. 

It is desirable to develop additional methods for studying processes with ab- 
sorbing states. In this work we study models that admit an active (nonabsorbing) 
stationary state in the infinite-size limit, but which always, for a finite system size, 
end up in the absorbing state. The quasi- stationary (QS) distribution for such a 
system provides a wealth of information about its behavior. (In fact, simulations of 
"stationary" properties of lattice models with an absorbing state actually study the 
quasi- stationary regime, given that the only true stationary state for a finite system 
is the absorbing one.) 

In particular, it would be valuable to have a simulation method that yields quasi- 
stationary properties directly. Currently available methods involve a somewhat com- 
plicated procedure for determining QS properties: a large sample of independent 
realizations are performed, and the mean <f>(t) of some property (for example the or- 
der parameter) is evaluated over the surviving realizations at time t. At short times 
times (f>(t) exhibits a transient as it relaxes toward the QS regime; at long times it 
fluctuates wildly as the surviving sample decays. Normally one is able to identify 
an intermediate regime free of transients and with limited fluctuations, which can 
be used to estimate the QS value of 0. This method requires careful scrutiny of the 
data and is not always free of ambiguity [11]. 

A number of strategies have been devised to circumvent the difficulties associated 
with simulating models exhibiting absorbing states. One involves replacing the 
absorbing state with a reflecting boundary in configuration space [12]. In another 
approach the activity is fixed at a nonzero value, in a "constant coverage" simulation 
[13] or a "conserved" version of the model [14]. If one includes a weak external 
source of activity, h, the zero-activity state is no longer absorbing, but it is possible 
to obtain information on critical behavior by analyzing scaling properties as h — > 
[11]. A further possibility is to clone surviving realizations of the process, enriching 
the sample to compensate for attrition as the survival probability decays [15]. While 
all of these approaches are useful, none affords direct, unbiased sampling of the QS 
state of the original process. 

For models without spatial structure, such as uniformly distributed populations 
or well stirred chemical reactors, the full QS distribution can be found from the 
master equation, via recurrence relations or an iterative scheme [16,17]. In models 
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with spatial structure, typified by nearest-neighbor interactions on a lattice, mean- 
field like approximations to the QS distribution have been derived, but descriptions 
in terms of one or a few random variables cannot capture critical fluctuations. The 
simulation method developed here does not suffer from this limitation. It provides 
a sampling of the QS probability distribution just as conventional Monte Carlo 
simulation does for the equilibrium distribution. 

We begin reviewing the definition of the quasi-stationary distribution. Consider 
a continuous-time Markov process X t taking values n — 0,1,2, S, with the state 
n = absorbing. We use p n (t) to denote the probability that X t = n, given some 
initial state X . The survival probability P s (t) = J2 n >iPn{t) is the probability 
that the process has not become trapped in the absorbing state up to time t. We 
suppose that as t — > oo the p n , normalized by the survival probability P s (t), attain 
a time-independent form. The quasi-stationary distribution p n is then defined via 

with p = 0. The QS distribution is normalized so: 

EPn=l- (2) 

n>l 



We now discuss the theoretical basis for our simulation method. The QS dis- 
tribution is the stationary solution to the following equation of motion [16] (for 
n > 0) 

dq n . . (0 s 

— = -w n q n + r n + r q n , (3) 



where w n = J2 m w m,n is the total rate of transitions out of state n, and r n = 
J2 m w n ^ m q rn is the flux of probability into this state. To see this, consider the master 
equation (Eq. (3) without the final term) in the QS regime. Substituting q n (t) = 
P s (t)p n , and noting that in the QS regime dP s /dt = — r = — P s J2m w o,mP m , we see 
that the r.h.s. of Eq. (3) is identically zero if q n = p n for n > 1. The final term in Eq. 
(3) represents a redistribution of the probability r (transfered to the absorbing state 
in the original master equation), to the nonabsorbing subspace. Each nonabsorbing 
state receives a share equal to its probability. 

Although Eq. (3) is not a master equation (it is nonlinear in the probabilities q n ), 
it does suggest a simulation scheme for sampling the QS distribution. In a Monte 
Carlo simulation one generates a set of realizations of a stochastic process. In what 
follows we call a simulation of the original process X t (possessing an absorbing 
state) a conventional simulation. Our goal is to define a related process, X£, whose 
stationary probability distribution is the quasi- stationary distribution of X t . (Note 
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that in order to have a nontrivial stationary distribution, X t * cannot possess an 
absorbing state.) The probability distribution of X t * is governed by Eq. (3), which 
implies that for n > (i.e., away from the absorbing state), the evolution of X t * is 
identical to that of X t . (Since Eq. (3) holds for n > 0, the process X t * must begin in 
a non-absorbing state.) When X t enters the absorbing state, however, X t * instead 
jumps to a nonabsorbing one, and then resumes its "usual" evolution (with the same 
transition probabilities as X t ), until such time as another visit to the absorbing state 
is imminent. 

A subtle aspect of Eq. (3) is that the distribution q n is used to determine the 
value of X t * when X t visits the absorbing state. Although one has no prior knowledge 
of q n (or its long-time limit, the QS distribution p n ), one can, in a simulation, use 
the history X* (0 < s < t) up to time t, to estimate the q n . This is done by saving 
(and periodically updating) a sample ni,ri2, ■■■,riM of the states visited. (If the 
state space is characterized by a small set of variables, one may instead accumulate 
a histogram H(n) giving the time spent in state n.) As the evolution progresses, X* 
will visit states according to the QS distribution. We therefore update the sample 
{rii, ri2, ...,um} by periodically replacing one of these configurations with the current 
one. In this way the distribution for the process X t * (and the sample drawn from 
it), will converge to the QS distribution (i.e., the stationary solution of Eq. (3)) 
at long times. Summarizing, the simulation process X t * has the same dynamics as 
X t , except that when a transition to the absorbing state is imminent, X t * is placed 
in a nonabsorbing state, selected at random from a sample over the history of the 
realization. In effect, the nonlinear term in Eq. (3) is represented as a memory in 
the simulation. 

To explain how our method works in practice, we detail its application to the 
contact process (CP) [4,5,8]. In the CP, each site i of a lattice is either occupied 
(cTj(t) = 1), or vacant (o"«(t) = 0). Transitions from <7j = 1 to <7j = occur at a rate 
of unity, independent of the neighboring sites. The reverse transition can only occur 
if at least one neighbor is occupied: the transition from o~i = to o~i = 1 occurs at 
rate Ar, where r is the fraction of nearest neighbors of site i that are occupied; thus 
the state o~i = for all i is absorbing. (A is a control parameter governing the rate 
of spread of activity.) The order parameter p is the fraction of occupied sites. 

To begin we study the contact process on a complete graph (CPCG), in which 
the rate for a vacant site to turn occupied is A times the fraction of all sites that 
are occupied, rather than the fraction of nearest neighbors. Since each site interacts 
equally with all others, all pairs of sites are neighbors, defining a complete graph. 
(The critical value in this case is A c = 1; the stationary value of p is zero for A < A c .) 
The state of the process is specified by a single variable n: the number of occupied 
sites. This is a one-step process with nonzero transition rates 
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on a graph of L sites. In Ref. [16] the exact QS distribution for the CPCG is obtained 
via a set of recurrence relations. 

We simulate the quasi-stationary state of the CPCG by realizing the process 
corresponding to the transition rates of Eqs. (4) and (5) and maintaining a list of 
M = 10 4 states. Each time a (nonabsorbing) state is visited, we update the list with 
probability 7 At where At = l/w n is the mean duration of this state. Whenever a 
transition to the absorbing state (n = 0) is generated, we instead select a state from 
the list. The results for a system of 100 sites, using 7 = 0.5, are shown in Fig. 1, 
illustrating excellent agreement with the exact QS distribution obtained in Ref. [16]. 
The simulation result is in good agreement with the exact result even for A = 0.1, 
deep in the subcritical regime, in which the lifetime of the original process is very 
short. 

Next we turn to the one-dimensional contact process (i.e, the model defined 
above on a ring of L sites). Although no exact results are available, the model has 
been studied intensively via series expansion and Monte Carlo simulation. In the 
QS simulations we use a list size M = 2 xlO 3 - 10 4 , depending on the lattice size. 
The process is simulated in runs of 10 s or more time steps. As is usual, annihilation 
events are chosen with probability 1/(1 + A) and creation with probability A/(l + A). 
A site is chosen from a list of currently occupied sites, and, in the case of annihilation, 
is vacated, while, for creation events, a nearest-neighbor site is selected at random 
and, if it is currently vacant, it becomes occupied. The time increment associated 
with each event is At = l/N occ , where N occ is the number of occupied sites just prior 
to the attempted transition [8]. 

In the initial phase of the evolution, the list of saved configurations is augmented 
whenever the time t increases by one, until a list of M configurations is accumulated. 
From then on, we update the list (replacing a randomly selected entry with the 
current configuration), with a certain probability p rep , whenever t advances by one 
unit. A given configuration therefore remains on the list for a mean time of M/p rep . 
(Values of p rep in the range 10~ 3 — 10~ 2 are used; the results appear to be insensitive 
to the precise choice.) 

Fig. 2 shows that our method reproduces the order parameter p obtained via 
conventional simulations. In Fig. 3 the QS and conventional results for the moment 
ratio m = {p 2 )/(p} 2 are compared. As discussed in Ref. [16], the lifetime of the 
QS state is given by r = l/p 1 . The lifetime obtained in QS simulations compares 
well with the lifetime obtained via conventional simulation (using P s ~ exp(— t/r)), 
as shown in Fig. 4. Detailed comparison shows that there is no significant differ- 
ence between the QS and conventional simulation results. For A < 3, conventional 
simulations are subject to relatively large uncertainties, since almost all realiza- 
tions become trapped in the absorbing state before the quasi-stationary regime is 
attained. At the critical point, quasi-stationary simulations require about an order 
of magnitude less cpu time than conventional simulations, to achieve results of the 
same precision. (Below A c the savings is even greater; well above A c the two methods 
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are equally efficient, since visits to the absorbing state are extremely rare.) 

As a further application of our method we study the critical CP (A c = 3.297848) 
on rings of 80, 160, 1280 sites. We perform between 5 and 30 realizations, each 
extending to 5 x 10 8 time steps. Results from the first 10 7 time steps (2 x 10 7 for 
L = 1280) are discarded from the averages; convergence to the QS state occurs 
on a considerably shorter time scale. Varying the list size M and the replacement 
probability p rep , we obtain results consistent with conventional simulations for M = 
1000 or greater, and p rep < 10 -2 . Smaller list sizes and/or larger values of p rep yield 
results slightly, but significantly, different from previous studies. On the other hand, 
there seems to be little point in using a list size much greater than the number of 
reinitializations. In practice the best approach is to use M = 10 3 or more, and to 
perform simulations with several values of p rep , to verify that the results show no 
significant dependence on this parameter. 

Our results permit us to refine the estimate for the moment ratio m = (p 2 )/{p) 2 
(p denotes the fraction of occupied sites). This quantity is analogous to Binder's 
reduced fourth cumulant [18] at an equilibrium critical point: the curves m(A, L) for 
various L cross near A c (the crossings approach A c as L increases), so that m assumes 
a universal value m c at the critical point. In Ref. [19] the estimate m c = 1.1736(2) 
for the critical one-dimensional CP was obtained from simulations of systems of up 
to 320 sites. Using the present method we reproduce and refine the results of [19]. 
Based on our data for L < 640 we estimate m c = 1.17370(2). 

At the critical point, the order parameter is expected to decay as a power law, 
p ~ L~ l3u± . The QS simulation data (for L = 80 - 1280) follow a power law 
with fiv± = 0.2525(5), compared with the literature value of 0.2521. The lifetime 
T ~ L v wl Vl - at the critical point; our simulation results yield v\\/i'± = 1.576(2) while 
the standard value is 1.5807. We find that including a correction to scaling scaling 
term of the form r ~ L u ^^ u± (l — bL~^) yields a substantially better fit (the sum of 
the squares of the errors is reduced by a factor of 2) if we use = 0.2. The best-fit 
parameters are b = 0.069 and v\\/v± = 1.581(2). (Precise determination of this 
correction to scaling term will require high precision data for a larger set of system 
sizes, a task we defer to future work.) Summarizing, the QS simulation method 
yields results fully consistent with conventional simulation, and with established 
scaling properties, when applied to the contact process on a ring. 

It is interesting to compare the QS distribution with that obtained using a reflect- 
ing boundary at n — 0, as was used in [12]. In the case of the CP on a complete graph 
the stationary distribution with a reflecting boundary (RB) is given in [16], where 
it is called the "pseudo-stationary" distribution. We find that for large systems, in 
the active phase (A well above A c ) the RB and QS distributions are essentially the 
same, but that nearer (and below) the transition the RB distribution yields a much 
higher probability for states near n = 1 than does the QS. The reflecting boundary 
is equivalent to a modified process X£ which, when a visit to the absorbing state is 
imminent, is always reset to the previous configuration. Since the latter is but one 
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step removed from the absorbing state, the buildup of probability in the vicinity is 
not surprising. In general, we should expect the QS and RB distributions to be in 
good accord when the lifetime of the process is reasonably long, since this implies a 
small QS probability in the vicinity of the absorbing state. 

In summary, we have devised and tested a simulation method for quasi-stationary 
properties of models with an absorbing state. The method is easy to implement, and 
yields reliable results in a fraction of the time required for conventional simulations 
in the critical regime, of prime interest in the context of scaling and universality. It 
also opens the possibility of investigating QS properties in the subcritical regime, 
which is essentially inaccessible to conventional simulations. We expect the method 
to be applicable to many problems currently under investigation, such as branching- 
annihilating random walks, conserved sandpiles, and stochastic population models. 
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FIGURE CAPTIONS 



FIG. I. QS distributions via recurrence relations (solid lines) and QS simulation 
(symbols) for the CP on a complete graph of 100 sites. A = 0.5, 1.0 and 1.5 (left to 
right). 

FIG. 2. Quasi- stationary density p in the one-dimensional CP. Open symbols: QS 
simulation, L = 20; filled symbols: QS simulation, L = 200. Solid lines represent 
results of conventional simulations. 

FIG. 3. Quasi-stationary moment ratio M in the one-dimensional CP. Symbols as 
in Fig. 2. 

FIG. 4. Quasi-stationary lifetime r in the one-dimensional CP. Symbols as in Fig. 
2. 
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